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We study double-barrier interfaces separating regions of asymptotically subsonic and supersonic 

■ flow of Bose condensed atoms. These setups contain at least one black hole sonic horizon from 
which the analog of Hawking radiation should be generated and emitted against the flow in the 
subsonic region. Multiple coherent scattering by the double-barrier structure strongly modulates 

■ the transmission probability of phonons, rendering it very sensitive to their frequency. As a result, 
resonant tunneling occurs with high probabihty within a few narrow frequency intervals. This gives 
rise to highly non-thermal spectra with sharp peaks. We find that these peaks are mostly associated 

QQ . to decaying resonances and only occasionally to dynamical instabilities. Even at achievable nonzero 

temperatures, the radiation peaks can be dominated by the spontaneous emission, i.e. enhanced 
zero-point fluctuations, and not, as often in analog models, by stimulated emission. 
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PACS numbers: 03.75.Kk Dynamic properties of condensates; collective and hydrodynamic excitations, 
superfluid flow 04.62.-|-v Quantum fields in curved spacetime 04.70.Dy Quantum aspects of black holes, 
evaporation, thermodynamics 
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Hawking radiation is one of the most intriguing yet unobserved predictions of modern physics. It is believed to 
^ ■ be generated near the horizon surrounding a black hole and to be responsible for its ultimate decay [H [1| . As noted 
by Unruh [3], a phenomenon analog to Hawking radiation should occur near a sonic horizon in a moving fluid, i.e. 
near the interface separating regions of subsonic and supersonic current. He argued that a thermal flux of phonons 
^ should be spontaneously generated from the sonic horizon towards the subsonic region. The effect finds its origin 
Q , in the impossibility of defining a global quasiparticle vacuum which suits both incoming and outgoing states. More 
O ' recently, it has been proposed that flowing condensates of bosonic atoms could provide interesting analogs of black 
hole physics [5|] and in particular of their Hawking radiation .6.- IO1] . The interface between subsonic and supersonic 
^ regions has also been shown to provide a scenario for the bosonic analog of Andreev reflection [ll| . 
^ \ Except for black hole lasers (discussed in Refs. d, [13, El, E^l ) ' most proposals predict a spectrum (phonon current 
' distribution per unit frequency) with a single peak at frequency a; = falling as l/w for small w > 0, where hw is the 
0^ ' quasiparticle energy measured with respect to the condensate chemical potential. This means that the low-frequency 
0^ , peak is essentially thermal in character. However, because of dispersive effects, the effective temperature characterizing 
' that zero-point radiation is no longer universal. In a Bose-Einstein condensate, the group velocity of phonons increases 
1 at high frequency, unlike in the original model considered by Unruh [l6j . As a result, this "superluminal" transport 
■ dilutes the sonic horizon into a spatial interval of flnite size [1, [l^ • The blue-shifting effect accompanying with the 
\ sonic horizon implies that the dispersive properties of phonons are always involved. Nevertheless, when the dispersive 
" . [ • length scale is much smaller than the horizon curvature scale, the temperature is determined by the local properties of 
condensate flow near the horizon (the gradient of the flow [3']) in strict analogy with the standard Hawking radiation 
which is fixed by the surface gravity of the black hole. This regime is found when the gradient is much smaller than one 
in units of the healing length Instead, when it is higher than one, the dispersion effects dominate and as a result 



d . the effective temperature is fixed by the healing length and the jump of the velocities across the sonic horizon [14 1. 
In any case, the effective temperature will be smaller than the chemical potential pL5| . As a consequence, a direct 
attempt to measure this radiation profile appears extremely difficult. An alternative proposal to indirectly measure 
Hawking radiation relies on the squeezed character of the state [1, 0] which could be observed at currently attainable 
temperatures by counting atoms on both sides of the horizon at coinciding times. However it should be pointed out 
that the main contribution to these correlations are due to the stimulated amplification of pre-existing phonons 0, Q , 
and not to the spontaneous amplification of vacuum fiuctuations. 

Here we study a new method which specifically aims at detecting the spontaneous contribution to Hawking radiation. 
This approach relies on the strong frequency-dependence of resonant tunneling through a double barrier structure. 
Such a sonic black-hole analog behaves as a Fabry-Perot resonator for quasiparticles, with the peculiar feature that 
quasiparticles propagate linearly against a condensate background which is itself governed by a nonlinear equation. 
The Hawking emitted phonon spectrum shows peaks at frequencies different from zero. We find that, at currently 
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achievable temperatures, thermal noise could be weak enough not to blur the characteristics of this resonant radiation, 
and a time-of-fiight experiment could allow for its detection. Our proposal is quite similar to the black-hole laser setup 
0, [13, [13, [13] in that sharp peaks are found in both cases. However, although dynamical instabilities may appear 
occasionally, they are not a necessary feature of these type of setups (see Section |V| . We propose here to focus on 
situations that are dynamically stable, i.e. where all peaks are due to resonances. The dynamical stability of the flow 
is likely to be a valuable asset in actual experiments. A systematic study of resonances and instabilities is currently 
under investigation. 

This paper is arranged as follows: In section II we present a mean-field study of the considered setup, discuss the 
main results and some general features. Section III is devoted to formulate the scattering problem of Bogoliubov 
quasiparticles propagating against the condensate background and to identify the essential features of Hawking radia- 
tion. In section IV we present and discuss the numerical results obtained for the Hawking radiation spectrum emitted 
from a double delta-barrier interface separating a subsonic and a supersonic region. Section V deals with the general 
distinction between quasinormal modes (or resonances) and dynamical instabilities, both of which are candidates to 
lie behind the sharp peaks in the radiation spectrum. Finally, Section VI is devoted to a summary and discussion of 
the main results. The main text is complemented by two appendices. Appendix A presents the analytical calculation 
of the mean-field model. Appendix B deals with the analytical resolution of the quasiparticle eigenvalue problem in 
the inhomogeneous region on the subsonic and supersonic sides near the double-barrier interface. 



II. FORMULATION OF THE MODEL. CONDENSATE WAVE FUNCTION. 



We study an atom transport setup which is schematically depicted in Fig. [T] A quasi-lD bosonic condensate, 
which occupies the left region j: < 0, is allowed to leak to the right though two identical delta potential barriers. The 
leftmost barrier is conventionally placed at a; = separated by a distance d from the second barrier. We will see that 
this double-barrier setup behaves as a resonant structure. For convenience, we neglect quantum fluctuations of the 
condensate stemming from its one-dimensional character [l3, [13] • 

We may decompose the Heisenberg second-quantized field operator 



$(a;, t) = e-'^*/''*o(a;) + S^{x, t) 



(1) 



into a stationary condensate wave-function e^*^*/''\['o(a;), with ^ the chemical potential, and its fluctuations 5'^{x,t). 
In this section we focus on the condensate behavior. At low temperatures and densities, the mean-fleld equation which 
governs the stationary flow of a Bose-Einstein condensate is the time-independent Gross-Pitaevskii (GP) equation for 
the condensate wave function: 



2m dx^ 



^^oix) = . 



(2) 



We wish to study the effect of a potential consisting of two delta barriers of equal strength, although for comparison 
we will occasionally consider the single barrier case (see Ref. [13]). Thus we will assume an external potential VcKt{x) 
which takes one of the two following forms: 



, . _ ( Vi{x) = hcuzS{x) , . 

^--i(^)-\V2{x) = hcuz[S{x)+S{x-d)] ' ^"^^ 

where Cu = \/ gunu/m, with n„ = lim3:^__oo ]^o(a;)]^, is the speed of sound on the subsonic, upstream (x < 0) side, 
and z the dimensionless strength of each barrier. The healing length on the asymptotic upstream side is S,u = h/mcu- 

In both the single- and double-barrier case, solutions can be found in which the condensate velocity is supersonic on 
the right of the second barrier, and subsonic from — oo to some point in the vicinity of the leftmost {x — 0) barrier (see 
Appendix for a general discussion of these solutions). Those flow profiles must have one or more horizons, defined 
as points where the local condensate velocity equals the local speed of sound. This subsonic-supersonic scenario is 
the same which has been shown to display the bosonic analog of Andreev reflection, where the supersonic side plays 
the role of the normal fluid [llj. Thus one may generally speak of Andreev-Hawking processes when dealing with 
scattering events undergone by elementary excitations of condensate flow through a subsonic-supersonic interface. 

Such sonic analogs relying on condensate flow have a superluminal dispersion relation at higher frequencies and 
their horizons do not imply a strict causal disconnection among regions, but they are sufficient to produce Hawking 
radiation analogs (see however Ref. [2]| for a discussion of a scenario where horizons would not be needed). For the 
case of a single delta barrier the only horizon lies on the near left of the barrier. The double barrier case is richer: there 
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appears one horizon in the vicinity of the x = barrier (which may he on either side) and possibly one or more horizon 
pairs. The leftmost horizon, and the only one in the single delta barrier case, can be viewed as a black hole (BH) 
analog, because there the flow goes from the subsonic to the supersonic side. The possible additional pairs of horizons 
appearing in a double barrier structure are the analogs of white-hole/black- hole pairs, where a white hole (WH) is the 
time-reversed version of a BH. White hole analogs seem to be extremely difficult to generate experimentally. There 
is a debate in the literature on whether white holes are stable at all, a question to which linear stability analysis has 
not yet provided an answer [2^ . 

To clearly identify what constitutes Hawking radiation, the experiment should ideally be done in such a way that 
the condensate wave-function is stationary and asymptotically flat: p'{±oo) = 0, where p{x) = |\E'o(a;)P [23|. In 
Appendixl^we show that, with those homogenous boundary conditions, there is one or more mean-field solutions for 
two delta potential barriers, and just one for a single delta barrier. A sketch of a typical density for two delta barriers 
is shown in Fig. [TJ 

In Fig. [2]wc plot the region of the {z, q) plane for which stationary solutions with homogeneous boundary conditions 
exist to the double-barrier problem. We recall that z is the dimensionless parameter characterizing the strength of 
the two identical delta barriers, while q is the condensate momentum on the subsonic (upstream) side. We note that, 
for a given value of q, there is a minimum and a maximum barrier strength between which a solution is guaranteed 
to exist for that particular value of q. This contrasts with the behavior of the single barrier case, for which only one 
solution exists for a given value of q. Interestingly, the resulting line z{q) line for the single-barrier case coincides with 
the upper boundary of the shaded region in Fig. [51 

Figure [3] shows a representative sample of solutions in the (d, z) plane, where d is the distance between barriers. 
Each curve represents a value of the condensate momentum q. Inspection of Fig. [3] reveals that, for a given q, there 
is a finite interval of allowed z values. That range has been shown in Fig. [21 Figure [31 reveals that, for given z and q, 
a multiplicity of d solutions exist. The two lowest d values can be characterized by two different values of p{0), which 
we call Pmax and /Omin- The upper d solutions come also in pairs and are regularly separated by a distance difference 
equal to the period of the nonlinear oscillations between the two barriers, which is the same for both pmax and Pmin- 

An interesting feature is that, at small barrier separations, a single q solution exists for a given d and z. For higher 
d (slightly above 2 for the sample of curves shown in Fig. [3]), two q solutions exist for given d and z. For still higher 
d, three q solutions exist for given d and z, and so forth. Since the barrier separation d and the barrier strength z are 
expected to be experimentally adjustable parameters, the clear trend is that, the larger the distance, the higher the 
number of allowed stationary solutions, a fact which could translate itself into instabilities. 



III. BOGOLIUBOV ANALYSIS 



A general and detailed introduction to the Hawking radiation physics in Bose-Einstein condensates, within a 
Bogoliubov-deGennes description can be found in Refs. @, 13 113, ■ this section we mainly introduce notation 
while refer the reader to those works for a more complete presentation. 

The quantum fluctuation part introduced in Eq. (|T]), 6'^{x, t), can be subject to a canonical transformation resulting 
in the expansion 



6^{x,t) = e-'^'^'^Yl [u^{x)e-'^^% + v*{x)e'^^'^l 



(4) 



where Ui{x) and v*{x) are the components of the wave function of the bosonic quasiparticle created by jj . These u, v 
components satisfy the bosonic Bogoliubov-deGennes (BdG) equations: 



u{x) 




v(x) 





H 



H giD*o(a;)2 
-5iD^S(^)' -H 





u{x) 




v{x) 



(5) 



2m dx-^ 



- p + Voxt{x) + 2giu\-9Qix)\' 



while the 7,7^ operators satisfy, for real w. 



7i,7- 



— / dx[u*{x)uj{x) — v*{x)vj{x)] = i^iSi 



(6) 



where the normalization vi can be set to ±1. 

In Eq. ([S]) (7iD = 21i}ao/{ma\) is the effective ID coupling constant, where ao is the s-wave scattering length and 
a_L the transverse confinement harmonic oscillator length, and Vc-x_t{x) is an externally imposed potential. We note 
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that the sum over i in the first equation may be interpreted as including both Ui ^ 0, with Vi > 0, or both Vi ^ 0, 
with uji > 0. 

The role of the horizon is appreciated when doing a WKB-type approximation for the excitations (see Refs. [23.I25I]). 
Close to the horizon are the turning points of the low energy classical trajectories where the WKB-solution is not 
appropriate. By matching the solutions on both sides and assuming some general scale separation (see also Ref. 
[2a|). one can show that in the case of just one horizon the profile of emitted Hawking radiation approaches l/cj 
at low frequencies. This important result guarantees the (approximately) thermal radiation profile. Thus at low 
frequencies, and in the absence of other scattering obstacles, the zero-point radiation has a 1/w behavior which makes 
it in principle difficult to distinguish from truly thermal behavior. 

Due to the presence of delta barriers in our chosen configuration, WKB-type approximations cannot be used 
uncritically. However, thanks to a theorem on dark-soliton perturbation theory (see Ref. [13]) we know that we have 
at our disposal a complete set of solutions on the left of the x ^ barrier. In the fiat supersonic region {x > d) the 
solutions are even simpler to work out. We refer the reader to Appendix IB] for both cases, a; < and x > d. The 
only non-analytically solvable problem lies between barriers (0 < x < d), but this is a finite region where numerical 
integration of the BdG equations ([5]) requires only a moderate computational effort. 

The next step is to identify the relevant scattering states. This has been done for Bose-Einstein condensates in 
studies of BH analogs and of bosonic Andreev reflection [ll|. In this article we focus on the scattering states and 

their connection to quasi-normal modes (or resonances) and dynamical instabilities. A systematic study of complex- 
energy eigenmodes is left for a future work [28J . In the asymptotic regions, propagating modes obey the Bogoliubov 
dispersion relation. Following Ref. we label the asymptotic regions with indeces u, for upstream {x — > — oo), and 
d for downstream (x — > c»). The upstream dispersion relation is 

LJu{k) = Vuk±Cu\kWl + ik^u)y^: (7) 

where Vu is the upstream flow velocity, and similarly for downstream ujd{k). A graph of this relation is shown in 
Fig. |3]for both the subsonic and supersonic side. The branches shown in blue/red correspond to the -f /— of Eq. ([7]) 
and can be shown to lead to positive/negative normalization. Modes are named after the sign of their group velocity 
(in/out according to whether they approach/leave the scattering structure), location (u/d) and, in the supersonic case, 
1-2 stands for modes with normal/anomalous (i.e. positive/negative) normalization. Importantly, the anomalous d2 
modes exist only for frequencies lo < Wmax- 

Here we are mainly interested in the scattering state with frequency w characterized by the incoming channel d2-in. 
Its wave function reads 



Ud2-in,a; 
1'd2-in,w 



(x) 
ix) 



-ipuix), X — > — OO 

ipd{x), X ^ oo 



(8) 



where 



i^d{x) 



■"d2(fcd2-in)e'«^" 
.l'd2(fcd2-in)e-*9<*^ 



»fcd2-i: 



v/27r|wd2(fcd2-in)| 



+ S'dld2('^) 
+ <S'd2d2('^) 



' Wdi(fcdi-out)e"?''^ ■ 
Wdi(fcdi-out)e~*«''^ 

' Wd2(fcd2-out)e"?''^ ■ 

Wd2(fcd2-out)e"'«''^ 



A/27r|K;di(fcdi-out)| 

^27r|wd2(fcd2-out)| ' 



(9) 



i>u{x) = 5ud2(w) 



^u(^u— out)^ 
Wu(Au-out)e"''?"^ 



pi/Cu _ out ^ 



\/27r|w„(fcu-out)| 



(10) 



where qu,qd are the up- and down-stream condensate momenta, w = duj/dk denote the relevant group velocities, and 
the ^'-matrix elements are shown. In Eqs. ([8l)- (|10|) Ua{k) and Va{k) are the spinor components for a quasiparticle of 
the uniform Bose gas, in channel a at momentum k. A similar decomposition can be made for the other in-modes. 
We note that the fc's of the previous equation are solutions of Eq. (O for a given uj lying in the interval < cj < ujmax, 
each solution defining a particular mode. Not shown in the previous equation but necessary for the matching, an 
evanescent solution exists on the subsonic side. The same holds on the supersonic side for lo > Wmax- To avoid double 
counting one can either choose all the positive normalization modes (for both in and out), as is done in Ref. [ll| . 
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and then one has to deal with negative frequency modes, or as is usual in this Hawking context, choose only positive 
frequency modes and then interchange 7^ O — 7I for the negative normalization ones (29j . In this work we adopt the 
convention of Wi > and write: 

/"OO 

S-^ix)^ du J2 i^iAx)li,^+vlJx)l\,J (11) 

Jo T ■ ^1 ■ 

7 ^u~in,al — 111 

(12) 

The normalization chosen in Eqs. [MTOl guarantees that modes are normalized to unit quasiparticle current and so 
[77, 7/' uj'] ~ ^ii'^i^ — ^')- An identical expression may be written changing in out. Standard scattering theory 
arguments show that the S'(aj)-matrix coefficients connect the in-modes to the out-modes: 

(7i-out,a;:7dl_out,td'7d2-out,c;) = (7u-i„,a; ' 7dl-in,tj ' 7d2-in,a;)<5^ (w) . (13) 

Due to the pseudo-Hermitian character of the BdG equations ([5]), this S'(aj)-matrix obeys, for frequencies w < Wmax, 
a pseudo-unitary relation S'' {uj)riS{ijj) — 77 with rj ~ diag(l, 1,-1). This enforces non-standard relations among the 
transmission and reflection coefficients, for example, |S'd2d2('^)P— |'S'ud2('^)P— |'S'did2(w)P = |<S'udi('^)P + |'S'didi(w)p — 
|5d2di(w)P = 1 (see Refs. giniEl). 

The case w > Wmax is simpler because the S{uj) matrix is unitary and there are no anomalous reflections or 
transmissions. Thus all asymptotic states can be chosen with v > 0. If we rewrite Ri{u!) :— S 1 1 (uj) , Tjj (uj) :— Sij{uj) 
with /, J different and taking values I,J = u, dl, then the usual unitarity relation \Ri{u})\'^ + |T/j(aj)p — 1 applies. 

The upstream phonon flux spectrum can be computed from these considerations and for to < Wmax shows the 
remarkable form [1, [13] : 

dIu~out I rt / .m2'^^u— in , i rr / _M2^-^dl — in , i r> / .m2 / ^-^u — ii 



^sMV=^ + \Sudii^r^^^^ + \s^M\'{=^ + i ■ (14) 

aco du! du! \ du J 

Assuming that we can populate the ingoing fluxes with comoving thermal populations, namely, dla^in/ du — 
riBi^a-in), whcrc a = u,dl,d2, 71.^(0) = (^e^^^ — is the Bose-Einstein occupation at the common tempera- 
ture /3^^ :— ksT, and fla-in is the comoving frequency of the mode a — in: 

f^a-in(t^) ^U - Vaka-iniuj), (15) 

where ka-iniju) is the solution to Eq. ^ for given lo, with Va = Vu,Vd the flow velocities. Equation (|14p reveals that, 
even at T = 0, a nonzero upstream flow of energy (phononic Hawking radiation) must be expected. The spectrum for 
cj > Wmax is of thc samc form as Eq. (jlip but with the last term removed, i.e. without any zero-point energy flux. 



IV. HAWKING RADIATION SPECTRA 



Figure |6] shows some frequency spectra of the upstream phononic flow [see Eq. (fT4)) ] for the setups and condensate 
solutions depicted previously in Fig. [5l with a one-to-one correspondence between the graphs in the two flgures. 
The top-left graph of Fig. [5] shows a structureless profile with a peak at w = followed by a l/ut tail, which is 
typical of Hawking radiation profiles in the absence of barriers, and very similar to what is obtained in the single 
delta barrier case. The message is that two nearby barriers behave similarly to a single barrier of double strength. 
We may also note that the zero temperature contribution (thick blue line) is not easily distinguishable from the total 
contribution at nonzero temperature (thin red line), because at low frequencies they both follow the same thermal 
law 1/uj. This property seems to be common to all structures which do not permit one or more nonlinear oscillations 
of the condensate between the barriers. An exception to this trend occurs when the delta barrier strength z is close 
to its upper limit. Then a bigger separation among barriers is possible and, as a consequence, a peak at nonzero cj 
may develop (see the upper-right graph at Fig. [B]) . A general trend that can be clearly appreciated in the rest of the 
graphs is that the largest the separation the more peaks appear. In particular the two bottom graphs in Fig. IHlexhibit 
a double peak in the allowed frequency interval. The reader might wonder why panels (e) and (f) of Fig. |6] look 
somewhat different since the parameters hardly vary, the only difference being a relative change in the inter-barrier 
distance of approximately 0.01. The reason is that the net amplification factor results from a rather complicated 
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expression involving interferences between the scattering coefScients on both sides of the scattering region. We refer 
to Eq. (69) and Fig. 5, right panel of Ref. 0, where a similar sensitive dependence has been found. 

We note that spontaneous phonons are generated at the event horizons and, after multiple scattering by barriers 
and horizons, they are partly emitted into the subsonic side. Those scattering events include normal processes {u and 
V components do not mix) and Andreev (or anomalous) processes {u converts into v or viceversa, see. Ref. [1 1] ) . The 
bigger the separation between the barriers, the more quasi-bound states can be accommodated between them and the 
more peaks appear in the Hawking radiation spectrum. 

An important point which will be discussed in greater depth in the next section, is that peaks can be due to reso- 
nances (also called quasinormal modes) or instabilities, so more information is needed to classify them. An unequivocal 
experimental signature of Hawking radiation (associated to zero-point quantum fluctuations of an otherwise stationary 
state) would undoubtedly be favored by a transport regime where the classical flow is dynamically stable. We recall 
that the S(ui) matrix is pseudo-unitary in the Hawking sector < cj < Wmax- The general spectral theorem on these 
matrices [30| guarantees that eigenvalues come in pairs of Si,l/s* or have unit modulus, \sj \ = 1. A general trend we 
have observed is that in the region w < Wmax/3 most of scattering matrices show only one unit-modulus eigenvalue. 
By contrast, in the region uj > Wmax/S, most of the iS-matrices have three unimodular eigenvalues. Nevertheless, the 
theorem mentioned before guarantees that the determinant will be a pure phase. Then, as in a conventional phase 
shift analysis of quantum mechanical scattering, a phase jumping upwards when crossing through a peak from low to 
high UJ can be interpreted as a resonance. On the contrary, a jump downwards will be reveal an instability. Hence, 
a most convenient way to detect resonances while distinguishing them from instabilities is to simultaneously plot the 
phase of the determinant of the S{uj) matrix. This corresponds to the thin green line appearing in all the graphs 
of Fig. [SI All curves show a clear resonance behavior, with the exception of the middle-right phase curve, which 
reveals an instability, as indicated by the sudden drop of the green line when traversing the instability. Most of the 
cases which we have studied show QNM behavior, but the the occurrence of instabilities is not so rare that it can be 
ignored. From our analysis, one cannot rule out the existence of strong instabilities with a large imaginary part in the 
eigenfrequency, since these would generate not easily detectable structures in the frequency spectrum. However, from 
the systematic analysis of Refs. where such instabilities where not found, one can conjecture that also in the 

present case those strong instabilities will not be found. An analysis of the time-dependent Gross-Pitaevskii equation 
in real time could establish that this is indeed the case. These considerations underline the need for a systematic 
search of instabilities and QNMs as poles of the propagator in the complex energy plane Ssj . 

The examples shown in Figs. [SJIHlwere chosen because they clearly reveal general trends. Unfortunately, in a time 
of flight (TOF) expansion, which correlates the long time density distribution with the momentum distribution of the 
initial state, their peaks barely stand out above the background of the depletion cloud, and even less when thermal 
fluctuations are taken into account. We have included in Fig. [7] a set of two more favorable setups which show large 
signals for zero-point Hawking radiation. The two bottom graphs show the momentum distribution of the initial 
trapped state and are computed in the approximation where only the subsonic flow is included and boundary effects 
are neglected. When a TOF measurement is performed in such systems, the depletion contribution is negligible at 
momentum values which however reveal clear resonant peaks. Zero-point Hawking radiation nearly exhausts those 
peaks at low temperatures and still gives the main contribution to the area under the peak at temperatures as high 
as 0.9/i. This fact could allow for the unequivocal detection of Hawking radiation. 

V. QUASI-NORMAL MODES (RESONANCES) AND DYNAMICAL INSTABILITIES 

A discussion of instabilities in BEG black- hole analogs can be found in Refs. [§, [13, [H, Hlj . For a study of QNMs 
in a similar context, we refer to Ref. [321. General studies of QNM in gravitational black holes and in optical analogs 
can be found in Refs. [3^ and [13, HH], respectively. A more complete study for the setup considered in this article 
will be presented [28| . 

We have said that peaks in the spectrum of Hawking radiation may be due to resonances or instabilities. The former 
are characterized by poles of the analytical continuation of the S{u}) matrix in the complex w-planc with Im(a;) < 
while the latter have Im(u;) > 0. Moreover, both types of complex modes may have |Im(a;)| so large that they are not 
clearly appreciated in the radiation spectrum, yet they can hide important instabilities. While a systematic search 
for poles is left for a future work, we discuss here some general features of the behavior of QNMs and instabilities in 
the complex /c-plane. 

When a linear stability analysis is made of a stationary solution of the GP equation instabilities appear as 
solutions of the BdG equations (O whose frequency uj has a positive imaginary part. Because complex frequency 
implies complex asymptotic wave numbers, these solutions must be localized in real space. For a given complex 
frequency, there are always four complex fc's (counting each different k with its multiplicity), two with Im(A:) > and 
two with Im(fc) < 0. 
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Like in the search for bound states in conventional quantum mechanics, a discrete, possibly empty set of modes is to 
be expected. The mode with the largest Im(a;) dominates the long time behavior and its inverse is a good measure of 
the decay time of the condensate due to that instability. There is a well defined procedure to quantize such unstable 
modes 36], which show up as essentially free particles (instead of oscillators). In open systems like the presently 
analyzed, there are however some subtleties in that quantization procedure which, to our knowledge, have not been 
properly addressed in the literature. Specifically, some spectral properties of the BdG operators are used which can 
only be guaranteed for finite-dimensional operators (37| . 

At a given real frequency between and Wmax we may have have four complex k solutions on each side of the 
interface. These are shown by circles in Fig. [H On the upstream side (left graph), they correspond to an incoming, 
outgoing, evanescent, and exploding solution. Downstream (right graph) we have two incoming and two outgoing 
solutions, corresponding to the normal (dl) and anomalous (d2) channels. For instance, a conventional retarded 
scattering state is characterized by one incoming channel matching all possible outgoing and evanescent solutions, 
with no amplitude for exploding or other incoming solutions. 

An instability is a bound state made exclusively of spatially decaying solutions that happen to match at a particular 
complex frequency (with positive imaginary part). Localized solutions correspond to Im(/c) < on the upstream side 
and to Im(fc) > on the downstream side. Figure [5] shows how the various k solutions evolve as uj varies from 
Ini(tj) = to Im(a;) > while Re(cLi) remains constant. Eventual instabilities can only be obtained by matching the 
evolved "u-out" and "evanescent" solutions upstream (Imfc < 0) and the evolved "dl-out" and "d2-out" downstream 
(Im/c > 0), with vanishing amplitudes for the other solutions. We emphasize that the wave function of these unstable 
modes involves exclusively analytical continuations of outgoing and evanescent solutions. 

QNM modes, on the other hand, are obtained from the analytical continuation to the Im(a;) < half-plane of the 
retarded Green's function of the time-dependent version of the Bogoliubov-deGennes equations ([S|) (see Refs. [13, 
for a discussion in an optical context but precisely translatable to the present one). More specifically, these QNMs can 
be obtained as an analytical continuation to Im(w) < of scattering states involving only outgoing and evanescent 
solutions, i.e. without the intervention of incoming waves. Like for instabilities, a discrete number of QNMs is to be 
expected. Therefore, in order to find resonances, we are only interested in the evolution of "u-out" and "evanescent" 
solutions on the upstream side, and "dl-out" and "d2-out" on the downstream side, ruling out the other solutions. 
Interestingly, this is exactly the same set of solutions which are relevant in the search for instabilities. We conclude 
that instabilities and resonances are obtained from the matching of the same set of solutions albeit in different regions 
of the complex energy plane. 

Hence, QNM/instabilities correspond to poles of the 5(0;) matrix when analytically continued to the lower/upper 
complex plane from the real energy line. This can be readily seen if one allows for a general coefficient in the in- wave 
component of the scattering solution Eq. for example. So both resonances and instabilities are candidates to 
explain peaks in the square of a given S{uj) matrix coefficient. The frequency of those possible underlying modes must 
of course have an imaginary part much smaller than their real part. Since the determinant of the S{lu) matrix is a pure 
phase (as noted in the previous section), it can be used to discriminate between both behaviors. Specifically, a jump 
upwards/downwards of the phase when traversing the peak with increasing frequency, implies a QNM/instability. 

Figure [SI together with other not shown data, reveals that peaks in the Hawking spectrum are mostly due to 
resonances, instabilities being more the exception than the norm. While this sampling is encouraging, a systematic 
search for QNMs and instabilities in various stationary fiows is left for future work. 

Finally we note that, both for QNMs and instabilities, the boundary conditions here described are different from 
those adopted in Ref . [3l|, . 

VI. DISCUSSION AND CONCLUSION 

We have studied the fiow of an atomic condensate through a double barrier interface separating regions of subsonic 
and supersonic flow. Such a setup provides a scenario where Hawking radiation into the subsonic side is enhanced 
at some frequencies due to multiple scattering of quasiparticles by the two barriers and the modulations of the 
condensate. The resulting highly non-thermal Hawking radiation presents peaks at frequencies which may lie well 
above the working temperature and thus can be unambiguously interpreted as stemming from quantum fluctuations 
of the quasiparticle vacuum. The non-thermal Hawking spectrum emitted by the double barrier interface represents 
and important advantage over the cases of single or zero barrier, where the low-frequency zero-point radiation has a 
thermal character which makes it more difficult to distinguish from genuinely thermal radiation. 

Our calculation has been based on a model of stationary fiow of condensate and quasiparticles through a double 
barrier structure with open boundary conditions, where the condensate density is asymptotically flat on both sides 
of the structure, while quasiparticle motion is described by scattering states characterized by incoming channels that 
are thermally populated. While a stationary scattering picture of transport has proved to hold predictive power in 
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electron systems, it still represents an idealized scenario in cold atom contexts, where stationary circuits have not yet 
been developed and where transport of finite-sized condensates is mostly investigated within a time-dependent scheme 
[38l |. As long as steady condensate transport is still an item for the future, it will be of interest to perform numerical 
simulations of time-dependent transport that may reveal those features of the Andreev-Hawking phenomena which 
we have explored here. 

An important question is whether, in currently achievable setups, the resonance frequency can be tuned to lie 
well above the currently attainable temperatures that would characterize the incoming quasiparticle population. We 
have seen that, while in a time-of- flight experiment the contribution from Hawking radiation peaks may be easily 
overshadowed by the contribution of the depletion cloud, setups can be designed where the resonant Hawking peaks 
are sufficiently sharp to be clearly visible even at temperatures comparable to the chemical potential. 

We thank A. Aspect, C. Diaz Guerra, L. Garay, P. Leboeuf, N. Pavloff, G. V. Shlyapnikov, and C. Westbrook 
for valuable discussions. This work has been supported by the joint France-Spain Accion Integrada HF2008-0088 
(PHC - Picasso Program). Support from MICINN (Spain) through grants FIS2007-65723 and FIS2010-21372, from 
Comunidad de Madrid through grant MICROSERES-CM (S2009/TIC-1476), and from the Swiss National Science 
Foundation, is also acknowledged. 



Appendix A: Mean-field calculation. 



In this and the next appendices we will use units such that ft = to = ^„ = 1, where ■= fi/ ^/TngiDnu is the 
asymptotic coherence length on the subsonic (upstream). This means in particular that the unit of frequency and 
energy is gmn-u- We will also rescale the wave functions by Then the equations in these units can be 

immediately read from Eqs. ©-([S]) by making h — m ~ guy ~ 1, which amounts to taking rt„ = lim^_j._oo |^o(2;)P = 1 
and Cu = I- The time- independent Gross-Pitaevskki equation © for the (scaled) condensate wave- function 'i'a{x) 
reads : 



*o(a:) =0 



(Al) 



We consider potential profiles Vcxt{x) consisting of one or two delta barriers, namely Vi{x) :— zd{x) and V2{x) :— 
z[6{x) + 6{x — d)]. The condensate will be asymptotically fiat and subsonic to the left of the first barrier and fiat and 



supersonic to the right of the last S barrier, 
reads [H: 



Then, we can choose phases such that for a; < the condensate profile 



^o{x) 
7(x) 



7(0) - iq 



\Jl — q^ tanh \/\ — q^{xQ — x) 



(A2) 



We note that xq ^ 0, which means that p'{Q ) < 0, where p{x) :— |'I'o(a;)P — > 1 as a: — ^ — cjo, see below. The chemical 
potential and uniform current read 



p = 1 



J 



q , 



(A3) 



where ^ q < 1 is required to have j > and subsonic fiow. 

On the supersonic side, the condensate wave function is a simple plane wave, which by current conservation has to 
be of the form: 



A ■ — 



lV8 



< 1 



(A4) 



where (j)a is a phase to be determined later. 
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There is a simple analogy which permits to qualitatively understand the behavior of the condensate wave-function. 
Shifting to an amplitude and phase representation, and splitting the GP equation (|A1I) into its real and imaginary 
components in the regions where Vc-xt{x) = 0, we may write 

= ^[A{xrcl^'{x)]^0. (A5) 

The third equation, which comes from the imaginary part of Eq. (jAll) . is but the continuity equation. The second 
equation, after multiplication by A'(x) and one integration, can be written as [cf. Eq. (jA3[) ]: 



E = -LL + v,{AM) 

Vi(A) :^ (l4)^' + |5-4;, (A6) 

E being an integration constant. Therefore (see e.g. Ref. [i^lSl) by looking at the position as a time coordinate 
and the amplitude as the position of a fictitious particle, this equation expresses the energy conservation for this 
particle under a force deriving from the potential Vq{A) (see left Fig. |9]). For q < 1 and A > this potential has 
a minimum at ^min [see Eq. ()A4p ]. on the left of a maximum at A^ax — 1- The flat solutions [A'{x) = 0] at 
the minimum/maximum are supersonic/subsonic [see Eq. (IA4[) and Ref. [l9j]. Integration of this equation with 
the initial condition lim2,_j._oo ^(a;) = 1 leads to (|A2[) . We note that the integration constant must be chosen as 
-Emax = Vq{Aynax) = l^g(l)- We wiU also usc the quantity £',nin = Vq(Aynin)- A delta barrier zS{x — xq) appears as 
an instantaneous kick, occurring at "time" xq, which adds a positive amount to the speed of the particle, which now 
becomes A'(xJ) = A'{xq) + 2zA{xo)- The "energy" of the particle changes instantaneously without changing the 
phase (j){xo) or "position" A{xo). In both the single- and double- barrier case, the first kick, S{x), must take place 
when the particle has negative "speed", ^'(0~) < 0, otherwise the kick would increase the energy with no chance of 
ending at ^min- The case of a single delta barrier at the origin, Vi{x), leads to a simple change in energy which can 
be solved to obtain z{q) [or q{z), as the dependence is monotonic, see Fig. 



-'-^max ^min — ^■^ ^min 



'^^^ = \^A^ ' 



min 



and the parameter xq can be extracted from (jA2[) by imposing 1^*0 (xg)! = A„ 

The case of two delta barriers, V2{x), is more involved. In this case, the energy after the first kick must lie between 
the two limiting values i?min < E < E^-^i^x, because if bigger than E^i^^, it would mean that the particle has acquired 
a positive speed, and the second kick could never send it to Ay^-m- To end up in A^-^ with zero kinetic energy, 
A'(d~) < is required. We conclude that, for a given current q, the maximum strength z which allows for a solution 
is (IA7[) . which is the unique z value for single delta-barrier. From the matching at a; = d it follows that the energy 
for the particle motion between kicks (0 < x < d) is 



+ Vq[A{x)] = En,in + 2z^A^i^ E{q, z) . (A8) 

Analysis of the x = Q kick leads to another equation which, combined with (jA8[) computed at a; = 0+, leads to the 
following two equations: 

E{q, z) - ^^^^ + Vq [Am = \+q' + 2zA(0)A'(0+) - 2^2^(0)^ (A9) 

which allows us to solve for A{Q)^A'{Q^) as a function of q, z. If we introduce the densities p{x) = A^{x),pQ ;= 
p(0),Pq := p'(0+), then Eq. (jA9l) can be rewritten: 
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po [E{q,z)-V,[^a)] 

+ 1/2 -E{q,z) 

2zpo , (AlO) 

z 

whose plot is shown in right Fig. [3] [45^. The minimum z which allows for a solution can be computed by enforcing 
that the two lines cross tangentially (see blue line in Fig. [2]). 

From the previous discussion, a method can be outlined to compute the possible distances d between the barriers 
for given q,z between the relevant limits: One takes one of the two solutions of (jAlOp . and propagates it (clockwise 
in the blue line of right Fig. [5]) until p{x) = ^^jn, p'ix) < 0. The "time" x needed is the minimum distance between 
barriers. An integer number of periods can be added, where the period is the time needed to get back to the initial 
point. Two families of solutions result, corresponding to the two solutions that can be chosen from Eq. (jAlOp . 

To obtain explicit formulae, it is necessary to solve for the motion between the delta barriers. A straightforward 
method is to integrate Eq. (|A8[) . written in the density representation: 




P'o ^ 



p'{x) 



2 



= (fix) - ly {p{x) - q^) - 2AE((7, z)p{x) 



4 

^E{q,z) := E^^^~E{q,z)>Q, (All) 

the last inequality requiring z < z(q), which will be assumed hereafter. The polynomial in p{x) on the r.h.s. of the 
first equation has three real roots ordered such that < ei < 62 < 1 < 63. The integral is tabulated (see Ref. [46l |) 
and the solution is: 



p{x) = ei + (62 - ei)sn^ {y/e^ - ax - Fq, k) 

F(arcsin(y^),fc), p',>Q 
2K(fc)-F(arcsin(yfer),fc), p'^ 



F, 



< 



V 63 - 61 



where F{(j),k) :— dO/ \Jl — fc^ sir? {9), K{k) :— F(7r/2,fc) are the incomplete and complete elliptic integrals of the 
first kind, and sn(M, fc) is a Jacobi elliptic function, whose period is 'iK{k). Some straightforward algebra allows to 
compute the possible inter-barrier distances in explicit form: 



a- Fo + 2nK(k) „ , , 

Ves - ei 
_ J ai, cn(al,fc)<0 
" \2K{k)-ai, cn(al,fc)>0 



ai 




^^sis — ^,fc|. (A12) 

62 - 61 



This identifies 2K{k)/ — 61 as the "period", i.e., the distance that can be added to the deltas to have another 
solution in the same family, as mentioned in the paragraph following the paragraph of Eq. (jAlOp . A plot of the 
possible d solutions for four possible values of the current is shown in Fig. [31 
Finally, the phase (f){x) can be computed from the last Eq. (|A5|) : 



(x) - m 



dx' 

q 



n am(Ve3 



61a; + Fq, fc), 1 , fc 

ei 



n ( am(Fo,A:),l - —,k 
ei 



(A13) 
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where, here and in the previous equations, cn(u, fc), am(M, k) are Jacobi elliptic functions, and n((/), m, k) := d9{{l~ 

TOsin^ 9)\/l — k"^ sin^ 9}~^ is the incomplete elliptic integral of the third kind. 
We summarize below the algorithms to compute the complete GP solution: 

Single delta barrier algorithm. Start with a given current < g < 1. Get the unique delta barrier strength, z, from 
Eq. (|T7| . Compute xo > from Eq. The solution is then Eq. (|X2|) for a; < and Eq. with (/)a = for 

X > 0. 

Double delta barrier algorithm. Start with a given current < q < 1. Compute the minimum delta barrier strength, 
z^in{q), as that which makes the two equations (IA10[) tangent in the po, p'o plane. Compute the maximum delta barrier 
strength, Zmax(<z) ~ z{q) from Eq. (|A7I) . Choose z^i-n{q) < z < z„^ax{q) and solve for po,p'o in Eq. (|A10p . Choose 
one of the two solutions. Find the roots of the polynomial Eq. (jAlip and order them ei < 62 < 63. The expression 
for p{x) between deltas can be obtained from Eq. (jA12p . Take one of the possible distances, d, from Eq. (jA12[) . 
Compute xq > from Eq. (|X2|) and the chosen value of pQ. The solution is Eq. (|X2|) for x < 0, Eqs. (|A12p . (|A13p 
with 0(0) for < x < d, and Eq. (|X4)) with (j)a = (l){d) from Eq. ((M3)) for x > d. 



Appendix B: Exact solution of BdG equations. 



1. Subsonic region 



For a general reference on the content of this Appendix, see Refs. [27l l47j . There are some sign conventions in 
those references which are not followed in the present work. It can be shown that a solution of the time-dependent 
GP equation (with no external potential). 



dt 2 dx 

makes the following system over-determined: 



d-^(x,t) ld^^(x,t) , , 

» ^^-^ ' +\^{x,t)\^'f{x,t) , (Bl) 
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wi{x,t) 






W2 (X,t) 




d 


wi{x,t) 




dt 


W2{x,t) 





-iX 'S*{x,t) 1 r ^1(0;,^) 
^'(a;,i) iX W2{x,t) 

i\^ix,t)\^/2 + iX^ -id^'^*{x,t)/2- X^*{x,t) 
i9^^'(x, t)/2 - X^ix, t) -i|*(x, t)| 2/2 - iX^ 



wi{x,t) 

W2{x,t) 



(B2) 



This is the so-called Lax pair of equations of the Zakharov-Shabat problem. The spinor solutions with certain 
boundary conditions of the first equation are denominated Jost functions. Here A is an auxiliary scattering parameter 
which allows to map the general solution of the Cauchy initial value problem posed by the GP equation into a linear 
scattering problem at t = from ^'(a;, 0). The time evolution of the scattering data is easily computed, and from the 
results at time t on can reconstruct 'i'{x,t) using inverse scattering techniques. Then, it can be easily shown that 
the squared Jost functions solve the linear perturbation problem associated to Eq. (IBip . which is described by the 
time-dependent BdG equations (see Refs. [39l - t4l| '): 



idt 



wj{x,t) 




_ wl{x,t) _ 





-a,,/2-f 2|*(x,i) 

~^*{x,t)^ 



^'(a;,t)2 
j2-2\^{x,t)\^) 



wj(x,t) 
wl{x,t) 



(B3) 



where we note the non-obvious order of the spinor components. 

In the particular case case of a stationary, asymptotically flat, and subsonic condensate solution, 'i>{x,t) 
e~'^'^*'^o{x), where ^o(a;) is given in Eq. (jA2p [see also Eq. (jA3p ]. a general stationary solution of (jB2l) is: 



wi {x,t) 
W2{x,t) 



^i{kx-et)/2 



g»(pt-gx-eo)/2^^^2;) 

e-'(''*-«^-^«)/2w2(x) 



(B4) 



where v^{—oo) = (i = 1, 2). The variables k/2,e/2 are introduced in such a way that k and e can be identified with 
the momentum and energy of the excitation after squaring. The phases are chosen to cancel the asymptotic phase of 
the condensate (except for a phase 1 — q'^ + iq) [see Eq. (|A2p ]. Substituting ansatz (IB4p into the first Eq. 
leads to: 
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d 


vi{x) 




dx 







7(a;)4 



-q)/2] j(x)- 
iq i[X — (fc - 



iq 



vi{x) 



(B5) 



which in the Umit x — > — oo and using the oo) = boundary condition, yields A: 



, q B 



(B6) 



Equation (jBSp can be solved by adding and subtracting both equations, which introduces S{x) :— vi{x) + 
V2{x),T{x) := vi(x) — V2{x). Solving for S{x) as a function of T(x),T'{x), and introducing the result into the 
other equation, leads to T"{x) + ikT'{x) ~ 0, whose only solution satisfying the boundary condition v'^{—oo) = 
[which implies r'(— oo) = 0] is a constant T{x) — Tq. Some further algebra leads to: 



vi{x) 



s{x) + T{x) ^ n 

2 2 

S{x) - Tjx) ^ To 
2 2 



«7(x) — k/2 



q/2 + X 
ij{x) — k/2 



q/2 + X 



(B7) 



To solve for the spectrum, e, we use the dt equation of (jB2l) . Due to the over-determinacy of Eqs. (jB2[) when t) 
is a solution of Eq. (jBll) . we only need one of them, e.g. the first one. Solving this equation in the limit x — t- — oo, 
the Bogoliubov spectrum is found to be [after invoking Eqs. (|B6p and (|A3[) ]: 



s = qk± k\ 1 + 



/c2 



where ± correlates with the same symbol in Eq. (|B6|) . From Eqs. (|B6 
that oo,t)|2 — |ti(— oo,i)|2 = 1, we can write the modes as: 



(B8) 

we get A + q/2 = e/fc, and demanding 



u{x, t) 
v{x, t) 



et) 



\e - kq) 



e 



2ij{x))) 
(l-±ik^2^J{x))y 



(B9) 



2. Supersonic region. 

In this case there is no need to use the technique of squaring the Jost functions, as the solution can be trivially 
worked out. We only quote the results, using the notation of Eq. (|A4[) and introducing p := q/Amm, c := Amin, which 
are, respectively, the condensate wave vector and the speed of sound in the supersonic region: 



u{x, t) 
v{x, t) 

u(x, t) 
v{x, t) 



^i{kx — et) 



_g-i(pa;-A't+0a)G_(fc^q) 

e-^ip^-i^i+'t>'^)G+{k,q) 



where G±{k,q) is a short-hand notation for 

G±{k,q) 



e/2 + _^ 1 

^fc2(fc2 +4c27 2 



1/2 



(BIO) 



(Bll) 



We are assuming, without loss of generality, that e > 0. The first solution ljBlOp is always valid and is normalized to 
|M(a;,t)|^ — |z;(a:,t)|^ = 1. The second solution is propagating only when e < cunuix and its normalization is anomalous: 
|u(a;,i)|^ — |f(a;,i)|'^ = —1. As expected, the spectrum is that of Bogoliubov quasiparticles for the supersonic region: 



e = pk ± 




(B12) 
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FIG. 1: A typical black-hole structure dis cusse d in this paper. The two vertical arrows represent the delta barriers. The profiles 
of the speed of sound [proportional to \J p(x), where p{x) is the local condensate density] and the speed of flow [proportional 
to 1/ p{x)\ are shown. Condensate flows to the right, while Hawking radiation is emitted into the left (subsonic) region. In 
this particular case, three event horizons exist since the two curves intersect at three points. The flow regime is such that the 
condensate density exhibits a depression between the barriers. 
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FIG. 2: For a structure of two delta barriers of identical strength, the region of z, q for which solutions exist is shown in grey 
color. The minimum and maximum delta barrier strengths z are plotted as a function of the current q (which in these units 
coincides with the upstream condensate momentum). Red curve shows the maximum z for a double delta barrier system, which 
is the same as the unique z value for a single delta barrier of the same strength [see Eq. (|A7I) ]. The blue line is the minimum 
z value allowing for a solution in the double barrier case [see Eq. (jAlOP and the ensuing paragraph]. 
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Z 

FIG. 3: Strength z of the double delta barrier as function of the inter-barrier distance d. The plotted lines ranging from 
right to left (with colors blue, red, green, brown) correspond to = 0.01, 0.05, 0.1, 0.3. Solid lines stand for solutions whose 
density profile undergoes zero or two oscillations between barriers, while the dotted lines correspond to solutions with one full 
oscillation. The multiplicity of the values of z as a function of d reflects the possibility of nonlinear oscillations between the 
barriers. The dots, with their respective labels, correspond to the various configurations shown in Fig. [6] (with velocity profiles 
given in Fig. [S]) and Fig. [T] 




FIG. 4: Dispersion relation [see Eq. ([7}] on the subsonic (left) and supersonic (right) sides. The blue/red branches correspond 
to positive/negative normalization as defined in Eq. ((6]). 
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FIG. 5: As in Fig. [T] plot of the local speed of sound (solid blue), proportional to p(x) and the local speed of flow (dashed red), 
proportional to Ij p(x), for a variety of setups. These correspond to relevant Hawking radiation profiles to be shown in the next 
Fig. El The parameters are for (a)-(f): number of added oscillations 0, 0, 1, 1, 2, 2; current = 0.05, 0.01, 0.05, 0.05, 0.05, 0.05; 
delta barrier strength z = 0.989, 5.851, 0.989, 2.012, 0.989, 0.989; of the two solutions for po in Eq. (|XTO)) . that with the smallest 
po is chosen in (a),(b),(e). The inter-barrier distance is d/^u = 0.445, 6.150, 2.871, 4.861, 5.161, 5.229. 
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FIG. 6: Plot of an upstream spectral radiation profile [see Eq. (|14p ]. so called Hawking radiation. The parameters of the various 
graphs are identical to their counterparts in the previous figure. Thick blue: current spectrum at ksT = 0. Thin red: finite 
temperature (in units of /i), 0.1 for the two uppermost curves and 0.3 for the rest). Dotted green: phase of the determinant of 
the S(ijj) matrix, hw-^^^j \i — 0.9299 for all the graphs except (b), which has &Jmax/pt = 0.9859. We draw the reader's attention 
to the fact that the peak in panel (d) corresponds to a dynamical instability, as can be seen from the drop of the green line as 
it approaches the peak for lower values of u, and not to a resonance, as it is the case in the other panels. 
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FIG. 7: Two setups of potential interest for a time-of-flight (TOF) experiment. Graphs in first/second row read as Figs. [5}l6] 
but with different parameters. The third row shows the momentum spectral resolution on the left (subsonic) region. The 
negative momenta correspond to atoms traveling towards the left. Sharp peaks are mostly due to resonant Hawking radiation. 
The background tail reflects atoms in the depletion cloud moving to the left. All graphs share the parameters q^^ = 0.1, 
hLOraax/ fJi. = 0.8612, oue added oscillation, and of the two solutions for po in Eq. (|A10[) . that with the largest value is chosen. 
They differ however in the barrier strength: z = 1.146, 0.968 for left and right column, respectively. In the second and third 
rows, the sequence of colors blue, brown, green, red, correspond to temperatures ksT/fj, — 0,0.3,0.6,0.9, respectively. The 
insets in the third row show the span of the peaks in greater detail. The determinant of the S{uj) matrix is not shown but 
reveals that the peaks are due to resonances. 
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FIG. 8: Migration in the complex fc-plane of the roots of the Bogohubov dispersion relation Eq. (O, for the two asymptotic 
sides (left is subsonic, righ is supersonic) of a BH configuration. Units of (left) or (right) are used. The paths shown 
result from tracking the k solutions as Im(a;) evolves from zero to a positive value, leaving the real part fixed. The labels 
and colors correspond to those of Fig. [4] but we have included the evanescent (green) and exploding (brown) solutions. The 
parameters are chosen so that the structure has a current q = 0.01 and Re(a;) is fixed to a value 0.8 < ojmax. 




>< 








^min(q) ^max(q)=1 

A 



^min(qr 



P(X) 



FIG. 9: Left graph: plot of the effective potential experienced by the fictitious particle presented in the discussion of Eq. (|A6|I . 
for a typical value of q between and 1. Right graph: plot of Eqs. (IA10|) (in blue the first, in red the second) for a typical case 
with an intermediate value of z. 



